clear all
figure
%Parameters%
rho=920; %kg/m^3
Q=50000; %J/mol
R=8.314; %J/mol/K
g=1.31; %m/s^2
%We now use a T dependent heat capacity
%We now use a T dependen thermal expansion coeff

%%%INPUTS%%%%%%
phi0=0.1; %initial porosity nominal 0.1
vrate=0.04; %nominal rate 0.04 m/yr (40 mm/yr)
%natron 1444 kg/m^3
rhosalt=1444;
fsalt=0.0; % nominal 0
ffsalt=0.0; %this is above ref
tcond=6000; % nominal 6 km
ncell=200;
etab=10^15; %Pa s %nominal 10 ^15
dip=pi*45/180; %ultimate dip in radians

Ts=100; %surface temp K (this is what Nimmo uses)
Tb=260; %basla temp K (this is what Nimmo uses)
rb=2*tcond; %bending radius
%%%%%%%%%%%%%%


dz=tcond/ncell;


%set up temperature profile
%and porosity and pressure
%This is set up so we have a perfect linear profile
T=zeros(1,ncell+2);
phi=zeros(1,ncell+2);
depth=zeros(1,ncell+2);
T(1)=100;
depth(1)=0;
depth(2)=dz/2;
for i=1:ncell
    if i>=2
    depth(i+1)= depth(i)+dz;
    end
    T(i+1)=Ts*(Tb/Ts)^(depth(i+1)/tcond);
end

depth(ncell+2)=tcond;
T(ncell+2)=260;

RR=rb+tcond/2-depth; % radial distance from bending center

P=depth*rho*g;
eta=etab*exp(Q/R*(1./T-1/Tb));
tau=eta./P;

%Porosity 65 Myr
phi=phi0*exp(-10^6*3.154*10^7*65./tau);


dt=3.154*10^7*1;
nstep=550000; 

phisv=zeros(ncell+2,nstep+1);
Tsv=zeros(ncell+2,nstep+1);
Dtop=zeros(1,nstep+1);
Dbot=zeros(1,nstep+1);
Dtop(1)=min(depth);
Dbot(1)=max(depth);
phisv(:,1)=phi;
Tsv(:,1)=T;


time(1)=0;
rhoav(1)=0;
%depth at time zero since diffusion is through slab
depth0=depth;



for j=1:nstep
    

    
    %update depth
    %remebetr dt in seconds vraste in m/yr
    theta=vrate*(time(j)+dt)/(3.154*10^7)/rb; 
    if theta<=dip
        depth=rb+tcond/2-RR*cos(theta);
    else
    depth=depth+vrate*dt*cos(dip)/(3.154*10^7);
    theta=dip;
    end
     %T dependent thermal conductivity 
     %This is from Petrenko and Whitworth 1999 (as refed by Nimmo and
     %Pappalardo 2016)
     %Cut off BCs for cell center temps and k
     for i=1:ncell+1
         k(i)=651*2/(T(i)+T(i+1));
         if i>=2
             %Temperature dependent heat capacity
             cc(i-1)=1925*T(i)/250;
         end
     end
     terms=1.0/rho/dz./cc; % for terms without k dT/dT = Delta q/ (c rho delta z)
                            % and Delta q = k Delta T/ delta z
     
     dTemp=diff(diff(T).*k./diff(depth0)).*terms*dt;
 for i=1:ncell
     %update temperature
     T(i+1)=T(i+1)+dTemp(i);
     
 end
    %update porosity
    phi=phi-phi.*P./eta*dt;
    phi(phi<0)=0;
    
    %update viscosity
    eta=etab*exp(Q/R*(1./T-1/Tb));
    
    %update pressure
    P=depth*rho*g;
    %update temp bc
    Treftop=min(260, Ts*(Tb/Ts)^(depth(1)/tcond));
    T(1)=(T(2)+Treftop)/2;
    
    
    %save things
    Tsv(:,(j+1))=T;
    phisv(:,(j+1))=phi;
    Dtop(j+1)=min(depth);
    Dbot(j+1)=max(depth);
    time(j+1)=time(j)+dt;
    
    %This is a reference temperature for same depth
    Tref=Ts*(Tb/Ts).^(depth/tcond);
    Tref(Tref>260)=260;
    %Reference viscosity
    etaref=etab*exp(Q/R*(1./Tref-1/Tb));
    tau=etaref./P;
    %reference porosity at 65 Myr
    phiref=phi0*exp(-10^6*3.154*10^7*65./tau);
 
    %only as much as thickness of conductive part
    % phi>phi ref negative anomaly 
    % T>Tref negative anomaly
    % alpha = 1.56E-4 T/250 (Kirk stevenson)
    alpha=1.56*10^(-4)*(T+Tref)/2/250;
    rhocell=0;
    for k=1:ncell
        n=k+1;
         
        %fraction of cell deeper than tcond
        %cells are at angle theta
         Dsalt=min(1,max(depth(n)+dz/2*cos(theta)-tcond,0)/(dz*cos(theta)));
         %fsalt-ffsalt ensures we dont double count ffsalt
         Dsalt=Dsalt*(fsalt-ffsalt)*(rhosalt-rho);
         
        rhocell=rhocell+rho*(Tref(n)-T(n))*alpha(n)+rho*(phiref(n)-phi(n))+Dsalt;
    end 
    rhoav(j+1)=rhocell/ncell+ffsalt*(rhosalt-rho);
    
    
end

%Plot if you want the nice maps
imagesc(phisv)
figure
imagesc(Tsv)

rhoav(1)=ffsalt*(rhosalt-rho);
figure 
 plot(time/(3.154*10^7)/1000,rhoav)


rhosum=0;
for i=1:numel(rhoav)
    rhosum=rhosum+rhoav(i);
    cumu(i)=rhosum/i;
end

figure
plot(time/(3.154*10^7)/1000,cumu)
